rm(list = ls())
library(Matching)
data("lalonde")
summary(lalonde)

# 估计倾向得分
glm1  <- glm(treat~age + I(age^2) + educ + I(educ^2) + black +
               hisp + married + nodegr + re74  + I(re74^2) + re75 + I(re75^2) +
               u74 + u75, family=binomial, data=lalonde)

# 为Match函数准备参数
Y  <- lalonde$re78 # 结果变量
Tr  <- lalonde$treat # 干预与否的标识

# 以倾向指数1对1的ATE匹配估计
X  <- glm1$fitted #倾向指数构造
rr  <- Match(Y=Y,Tr=Tr,X=X,M=1,estimand = 'ATE')
summary(rr)
# 检查数据平衡性
mb  <- MatchBalance(treat~age + I(age^2) + educ + I(educ^2) + black +
                      hisp + married + nodegr + re74  + I(re74^2) + re75 + I(re75^2) +
                      u74 + u75, data=lalonde, match.out=rr, nboots=500)
# 检查结果不错，t值较大
# std mean diff 即为标准化均值差异
# var ratio 即为二阶矩检验，未取对数而已
# KS检验原假设：两组样本来自的总体分布无显著差异

# 检查条件独立性假设
## 伪结果变量re75
rr  <- Match(Y=lalonde$re75,Tr=Tr,X=X,M=1,estimand = 'ATE')
summary(rr)
## 伪干预方法
if (FALSE){
  Tr <- lalonde$treat[lalonde$treat==0]
  Tr[sample(1:260,120)] <- 1
  rr <- Match(Y=lalonde$re75[lalonde$treat==0],Tr=Tr,X=X[lalonde$treat==0],M=1,estimand = 'ATE')
  summary(rr)
}
 


# 以马氏距离匹配1对1匹配
X <- as.matrix(lalonde[,1:8])
rr  <- Match(Y=Y,Tr=Tr,X=X,M=1,estimand = 'ATT',Weight = 2)
summary(rr)

# 若觉re75收入匹配效应不好，则可偏误纠正
Z <- as.matrix(lalonde$re75)
rr  <- Match(Y=Y,Tr=Tr,X=X,M=1,estimand = 'ATT',Weight = 2,Z=Z,BiasAdjust = T)
summary(rr)
# 检查数据平衡性
mb  <- MatchBalance(treat~age + I(age^2) + educ + I(educ^2) + black +
                      hisp + married + nodegr + re74  + I(re74^2) + re75 + I(re75^2) +
                      u74 + u75, data=lalonde, match.out=rr, nboots=10)
# 它的检查结果就不是很好，因为t值都是接近0的。(nboots=10)
